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Abstract. We present a variational algorithm for solving the classical inverse Sturm-Liouville problem 
in one dimension when two spectra are given. All critical points of the least squares functional are at 
global minima, which suggests minimization by a (conjugate) gradient descent algorithm. Numerical 
examples show that the resulting algorithm works quite reliable without tuning for particular cases, even 
in the presence of noise. With the right choice of parameters the minimization procedure converges very 
quickly. 



1. Introduction 

The one dimensional inverse Sturm-Liouville problem and methods for its numerical solution are very 
well studied [8, 11, 2]. Still there is no general purpose algorithm, which efficiently satisfies all needs (see 
eg. [12]). . ^ 

We present a variational algorithm which proves to be very robust under noisy input, does not 
require any special tuning or additional input besides the partial spectra, and is quite efficient in the 
class of variational methods. In section 2 we will define the functional and prove the absence of true local 
minima, which could trap our minimization procedure. Numerical examples, the optimal choice of the 
weights and the performance under noisy input will be discussed in section 3. Section 4 is concerned with 
implementation details and compares this algorithm to related approaches, and section 5 finally contains 
the proof of the linear independence of the squares of eigenfunctions. 

2. Definition and Properties of the Functional 

We consider the Sturm-Liouville equation 

— u" + q(x)u = Xu (SL) 
on [0, 1] with q(x) G £ 2 ([0, 1],K) real, and separated boundary conditions 

u(0) cos a + u'(0) sin a = 0, u(l) cos/3 + u'(l) sin/3 = . (a/3) 

For the asymptotics of the eigenvalues (A„)^L wrt. boundary conditions (a/3) there are three cases (see 
eg. [13, 4, 3, 9]): 

f lV -lsS + /o« ds + (I » ifsmasin/3^0 

An = < 7r 2 (n + 1/2) 2 + 2 Z S ( p-*f + Jo 1 ds + a " ifsmasin/3 = 0Asin( y 9-a) ^0 (2.3) 

[ 7r 2 (n + l) 2 + f Q q ds + a n if sin a = sin (3 — 
where («„) G I 2 . 
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It is a classical result, that two spectra determine the potential uniquely. Fix three angles a, 0, 
and 7 with sin(/3 — 7) 7^ 0. Let A g ^„, i e {1,2} denote the eigenvalues with respect to the boundary 
conditions (a/3) resp. (oij). 

Theorem 2.1 (Borg [1], Levinson [5]). Given potentials Q,qE i 1 [0, 1] with 
^Q,i,n = \,i,n for all n e No, i = 1,2, 

i/ien Q = q. 

Suppose we are given (partial) spectral data \Q,i, n with (i,n) in 7 C {1,2} x No of an unknown 
potential Q. For a trial potential q and positive weights cj^„, we define the functional 

G(q) = ^i,n{\,i,n - ^Q^nf ■ (2.4) 

(i,n)el 

If I is infinite and the sequence (w», n ) is summable, this converges because of the asymptotics of the 
eigenvalues (2.3) 

G(q) = V U i>n ( / Q- q + Ci >n ) < 00 , 
(i,n)€/ U ° ^ 7 

where Ci jn 6 i 2 ,i = 1, 2. 

If moreover 7 = {1,2} x No, the functional is zero if and only if q = Q (by theorem 2.1). In the case 
of partial knowledge of the spectrum, a solution of G(q) — includes all information given about the 
unknown potential Q. But we only have a chance to get close to the original potential if there is not too 
much information in the higher eigenvalues, i.e. if the potential to be recovered is sufficiently smooth [7]. 

An important question is now, for which sequences (Aj ; „) there is a q with G(q) — 0, i.e. Xq^ jn — ^i,n- 
It is easy to see with the classical methods [13], that it is necessary for the eigenvalues to interlace 

Al, n < M,n < M.n+l Or A 2! „ < \l t n < A 2 ,„+i . 

It is also known that if we choose either 

(i) sin a, sin 0, sin 7 ^ (Levitan [6]) or 

(ii) a = = 0, sin 7 ^ (Dahlberg Trubowitz [3]), 

then the interlacing property in connection with the correct asymptotics is sufficient for the existence of 
a potential q with G(q) = 0. 

An interesting feature of this functional is that all its critical points are at global minima. To prove 
this, we first compute the derivative of G. The derivative of A g; i in wrt. q in direction h is 



\,iA h \ = [ h 9q,i,n dx - 
Jo 



where g q ,i, n denotes the normalized eigenfunction corresponding to \ q ,i, n (sec [9] for a nice proof). Thus 
the derivative of G is given by 



G[h](q) = 2 V / Ui, n (\ qiit n - \Q,i tn )g 2 qin h&x . 

12 „^ r J0 



(i,ra)€l* 

We note that if ncji^ is summable, then the gradient 

VG(q) = 2 (»+lK»(\i,»-Ag,i,n)4T4 ( 2 - 5 ) 

(i,n)£l 

is in H 1 because ||5g,i,„||fl-i = 0(n) [13] and A ?>ii „ - A Qil! „ = 0(1). 

Theorem 5.2 below shows that the cigenfunctions g\i n are linearly independent in H 1 . This 
immediately implies the essential convexity of the functional: 
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Theorem 2.2. If I is finite or (nwj ;n ) is summable, the functional G has no local minima at q with 
G(q) > 0, i.e. 

G[h] (q) = V/i G H 1 [0, 1] <^=> G(q) = . 

Thus a conjugate gradient algorithm will not get trapped in local minimas, as we will also observe 
in the examples. 

Finally, we want to address the question, which of the infinitely many q with G(q) = our algorithm 
will select. Let us first look at a related functional, which actually was also our first try. It is given by 

= E ( X i-^ - X Q-"> + I - «) dx ) ' ( 2 - 6 ) 
(i,n)ei V Jo J 

with derivative 

d[h](q) = 2 V / f A,, i>n - A Qlj , n / (Q-^dz) l)/ida;. 

In case the mean of Q is known, this functional works as well as the other. Its gradient flow leaves the 
mean of q constant by construction. The function j| ( — 1 is the gradient of \ q ,i. n — J qdx and the 

direction of strongest increase of \ q ^. n which leaves J Q qdx fixed. 

Choosing a = (5 = and 7 = it/2, the asymptotics of the squared eigenfunctions are given by 

2 = ( 1 - cos ((2n + 2)ttx) + Oin- 1 ) ,i=l 
9 i' l - n \ 1 -cos((2n+ l)irx) +0(n- 1 ) ,i = 2 " 1 j 

Hence the functions {g^ i n — l|n € No,i = 1,2} U {1} arc almost orthogonal for large n. Now, the 
derivative of our functional (2.4) can be written as 

G[h]{q) = 2 f Wj,„(A g ,i,„- \Q,i, n ){{gli, n -l) + l)hdx, 

and the difference \ q + c ,i.n — f(q + c)dx = \ q ,i. n — J qdx is invariant under adding a constant function 
c to q. It follows that a gradient flow of our functional G leaves \ q ,i, n — Jq qdx with (i, n) ^ I almost 
invariant for n large enough. 

Since also a conjugate gradient descent is just an approximation of the gradient flow, in practice 
we do get some little changes in the higher eigenvalues. A similar argument holds for the case 
sin a, sin [3, sin 7 7^ 0, but there the asymptotics for i = 1,2 are equal. So we still have asymptotic 
orthogonality in n, but can not separate i = 1 and i = 2. 



3. Numerical Examples 

We use the standard Polak-Ribiere conjugate gradient descent algorithm [10] to approximate the gradient 
flow and thus minimize the functional. To give the basic idea, we explain the simpler steepest descent: 

(i) choose initial potential q = q 

(ii) while G{qi) too big do 

(a) compute the gradient VG(qi) 

(b) minimize the one dimensional function G(qi — aVG(qi)) wrt. a 

(c) set qi + i equal to the minimizing potential 
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Figure 1. Plot using u){ n — 1. The right two graphs show G{(^) 
and A2 over the number of iterations. 




Figure 2. Plot using uji t „ = (n + 1) 2 . On the right side of the A2 plot we zoom in to show the 
convergence process in more detail. 

This straight forward minimization scheme has a major disadvantage: consecutive gradients are always 
orthogonal. To avoid this, conjugate gradient descent computes the direction for the one dimensional 
minimization using the current and previous gradients. 

The first natural question is the optimal choice of the constants u>i_ n . A short look at the leading 
coefficients of the asymptotics (2.3) may suggest something like Wj >n = (n + 1)~ 2 . But the leading term 
of \,i, n — ^Q,i,n is f (q — Q) dx = 0(1), which corresponds to choosing u) itn = 1 and yields the fastest 
convergence in all examples. 

Figures f and 2 were computed with Wi >n — 1 resp. u>i jrl = (n + 1)~ 2 . For all plots (except 2 and 3) 
we chose the lowest number of iterations after which the plot stabilizes, while usually already the second 
iteration reveals the global structure of the potential. By default we use 30 pairs of eigenvalues, qo = 0, 
Wi tn = 1, and the boundary conditions a = (3 = 0, 7 = it/2. We give the current value of the functional 
G(q) as well as the L 2 error A 2 = ||g — QW2, and the maximum of the difference of the eigenvalues 
A\ = max( i; „\ e/ {|A 9i i ; „ — Aq^ : „ |}. In addition if there is a visible difference, we also plot the original 
potential Q using a dashed line. 

We see that both G(q) and A2 converge much faster for u>i tn = 1. With a>j )Tl = (n + 1)~ 2 on the 
other hand, q moves through smoother functions (as in figure 2). In contrast to the optical impression 
these are in general worse approximations wrt. G(q) and A 2 . But finally q, after around 100 iterations, 
will also converge to the function shown in figure 1. From the good average convergence of G(q) in both 
cases, we can tell, that our algorithm also in practice does not get near a local minimum, as proved in 
theorem 2.2. 

The optimal choice of boundary conditions can already be guessed from the asymptotics of the 
squared eigenfunctions (2.7). Choosing a — (3 — and 7 = ir/2, the functions g^i^ — 1 are almost 
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Figure 4. Plots with added white noise of absolute value smaller 
or equal 0.01 resp. 0.1. 



orthogonal for large n and i = 1,2 (c.f. (2.7)). In contrast, if sin a, sin/3, and sin7 are all non zero, 
9q,i,n ~ 1 an d 9q,2,n ~ 1 wm S e t close for large n. Therefore, in the latter case, the algorithm will converge 
much slower. This can also clearly be observed in numerical examples, like in figure 3. But also there, 
after approximately 620 iterations, we finally will get a solution which is close to figure 1. 

Another important aspect for applications is stability against noise in the given spectral data. In 
figure 4 we computed two examples with random noise (Aq^^ — Xq^^I < r, with r = 0.01 and r = 0.1, 
respectively. It is remarkable, that the convergence speed in G(q) is not significantly affected by the 
random noise; both examples reach G(q) ~ 10~ 18 in 30 iterations. 

For testing the robustness, we also fed the functional with the following random sequence 

9.99742, 11.6265, 14.4527, 23.9247, 26.2413, 31.091, 40.6658, 48.1088, 53.5093, 60.9088, 

which is far from the generic asymptotics. Convergence of G(q) is slow but steady. In log scale the graph 
of G(q) over the number of iteration looks similar to those given above. From G(q) — 36287 initially, we 
get down to G(q) = 3.56528 • 10~ 9 in 450 iterations. 

Finally figure 5 shows the results for the other examples of [2, 11]. 
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Figure 5. Continuous examples. 



4. Implementation and Comparison 

The implementation of the conjugate gradient minimization algorithm is straight forward (see eg. [10]), 
but our application crucially depends on a reliable eigenvalue problem solver. After some trial and error 
we settled for d02kdf from the Numerical Algorithms Group (NAG). For standard numerical routines we 
used the GNU Scientific Library. 

Another important point is the internal representation of the functions. We used cubic splines on 
2000 intervals of equal length. For 30 pairs of eigenvalues we already get satisfactory results at 180 
intervals. But since the computation time depends only roughly linear on this number, we did not study 
its influence on the quality of the approximation. 

Because all high level routines, like the eigenvalue solver, only sample these functions on some points, 
the details of the representation are not that important. But since very many values are needed, this is 
the part in the algorithm, where speed optimizations are most valuable. 

Compared to the variational algorithm by Brown et al. [2], which uses other spectral data, our 
functional is more expensive to compute, because they only have to solve initial value problems, whereas 
the present algorithm requires solution of the eigenvalue problems. On the other hand they often need 
about 10-100 times as many iterations to get similar results. Another plus on our side is, that even in 
presence of noise, we know to get closer to our goal each step and do not have to regularize the process 
by limiting the number of iterations. 

The method of Rundell and Sacks [11] uses the same spectral data as our algorithm and is without 
any doubt much faster, but needs the mean Jq Q dx as additional input, which has to be guessed from 
the (partial) spectral data. Our method does this automatically. In the case of figure 1 for example, the 
error of the mean is 2 • 10~ 3 % and using 5 pairs of eigenvalues we still only get an error of about 0.2%. 
This is probably better than an independent algorithm for computing the mean value from spectral data 
could do. 

5. Linear Independence of Squared Eigenfunctions 

Define the Wronskian [/, g\ = fg' — f'g and the bilinear form 
r : H l {% 1],M) 2 — > R 

(/>S0 ^ Jol/'^da; 

which is bounded by 

|r(/,ff)|<||/lk 1 Nk,i.e.||r(/,.)|| = ||/lk 1 . 

(We use the definition = \J \\.f\\ 2 L 2 + H/'H^ with distributional derivatives.) In particular T is 

continuous on H 1 . We have the following rules for the Wronskian: 
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(i) If 9, hj] = gj[f, h] + fh[g, j] = fj[g, h] + gh[f, j] 

(ii) If the functions /i and f 2 fulfill the condition 

fi(a) cos a + /■ (a) sin a = i = 1,2, 
then [/i,/ 2 ](a) = 0. 

(iii) For two arbitrary solutions /i and / 2 of the equation (SL) with eigenvalue parameters Ai and A2 we 
have 

[a, ./y' = /1/2' - .f2.f1 = (Ai - A2X/1/2 . 

In this section we are only talking about a single g and therefore drop it from the index. Let s ijM and 
Ci tn be the solutions of the differential equation (SL) for the eigenvalue parameter A^„ and initial values 

Si,„(l) = sin/3, Ci,„(l) = sin 7, 



s: 



< >n (l) = -cos/3, c' in (l) = -cos 7 . 



It is well known and easy to see that the Wronskian of these solutions is constant and we can compute 
its value at 1: 

[si, n , Ci_r\ = - sin/3cos7 + cos/3sin7 = sin(7 - (3). 

The normalized eigenfunctions are g\, n = Si,n/||si,n||2 and g 2 . n = c 2in /||c 2)n || 2 . Now we prove the central 
lemma, which builds on the ideas of a similar result for the Dirichlet case in the book of Poschel and 
Trubowitz[9]. A result in the same spirit can also already be found in the paper of Borg [1]. 

Lemma 5.1. Given three angles a, (3, and 7 with sin(/3 — 7) 7^ and denote the L? -normalized 
eigenfunctions of the Sturm- Liouville equation with boundary conditions (a/3) and (cry) by gi ;n , i = 1,2. 
With Si^ n , Ci^ n as defined above we have the following relations for the squared eigenfunctions for all 
i,j E {1,2} and m,n E No-' 

ft) T (gln,gl m ) = 

(ii) T(c l . n s l . ni g 2 3m ) = (-If sin(7 - f3)5 n>m 6ij 
Proof, i) 



T (9i, n ,9i, m ) = 2 / 9i,ngi,m[9i,n, 9i,m] dx 

Jo 



If n = to, the term clearly vanishes. If n ^ to, we get 

2 f 1 

^i.9i : n-> 9i,rn) = 7 7 / [9i,n •> 9i,m] [9i,n •> 9i,m] dx = — - [9i,m9i,m] |q = 0. 

because g\ n and gi m satisfy the same boundary conditions (rule ii). 
ii) 



2 f 1 

JO 

- n, the first term vanishes anc 

/ si,mgi,m[ci, m , gi, m ] = I < 
Jo Jo 



If i = j = 1, to = n, the first term vanishes and we are left with 

9i, m [ci,m, si. m } = - sin(7 - (3) , 
and if i = j = 2, to = n, the second term vanishes and we get 

/ C2. m g2, m [ s 2,m,g2,m] = / #2 m [S2,m , C 2 ,m] =+ SUl(7 - (3) . 

Jo Jo 
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If (i, n) ^ (j, m) we compute 

2 1 f 1 

/ n,n — ^j.m Jo 



[ s i,ni 9j,rn\ [Q,n? 9j,rn] |q 



We note that A« ?n — \j im ^ in this case, because the eigenf unctions satisfy the same boundary conditions 
at and different boundary conditions at 1 . □ 

Theorem 5.2. With the notations of the above lemma, the set of squared eigenf unctions 

{ff? n |(*,n)e{l,2}xN } 
is linearly independent in H 1 . 

Proof. Suppose for some fixed (i,n) we have 

9i,n = }^ akf 3k 

in H 1 , where a k € R and g k = 9j k ,m k with {jk,m k ) ^ (i,n). But this would imply 

(-I) 1 sin(7 -13)= T(a tn s i>n , g? >n ) = T I c hn s^ n , ^ a k g\ J = ^ r(c i;n s i;n , a fc ^) = 

V fceNo / fceNo 



□ 



Acknowledgments 



The author wishes to thank Ian Knowles and Rudi Weikard for drawing his attention to this topic and 
for the interesting discussions during his time at the University of Alabama at Birmingham. 



References 

[l 
[2 



[3 

K 

[5 
[7] 
[8 
[9 
[10: 
[11 
[12 
[13; 



G. Borg. Eine Umkehrung der Sturm-Liouvillcschcn Eigcnwertaufgabe. Bcstimmung dcr Diffcrcntialglcichung durch 

die Eigcnwerte. Acta Math., 78:1-96, 1946. 
B. M. Brown, V. S. Samko, I. W. Knowles, and M. Marietta. Inverse spectral problem for the Sturm-Liouville equation. 

Inverse Problems, 19(l):235-252, 2003. 
B. E. J. Dahlberg and E. Trubowitz. The inverse Sturm-Liouville problem. III. Comm. Pure Appl. Math., 37(2):255- 
267, 1984. 

E. L. Isaacson, H. P. McKean, and E. Trubowitz. The inverse Sturm-Liouville problem. II. Comm. Pure Appl. Math., 
37(1):1-11, 1984. 

N. Levinson. The inverse Sturm-Liouville problem. Mat. Tidsskr. B., 1949:25-30, 1949. 

B. M. Lcvitan. Inverse Sturm-Liouville problems. VSP, Zeist, 1987. Translated from the Russian by O. Efimov. 
M. Marietta and R. Weikard. Weak stability for an inverse Sturm-Liouville problem with finite spectral data and 

complex potential. Inverse Problems, 21 (4): 1275— 1290, 2005. 
J. R. McLaughlin. Analytical methods for recovering coefficients in differential equations from spectral data. SIAM 

Rev., 28(l):53-72, 1986. 

J. Poschel and E. Trubowitz. Inverse spectral theory, volume 130 of Pure and Applied Mathematics. Academic Press 
Inc., Boston, MA, 1987. 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes in C. Cambridge University 

Press, Cambridge, second edition, 1992. 
W. Rundcll and P. E. Sacks. Reconstruction techniques for classical inverse Sturm-Liouville problems. Math. Comp., 

58(197):161-183, 1992. 

P. Schapotschnikow. Inverse Sturm-Liouville problems on trees, a variational approach. Master's thesis, School of 

Computer Science, Cardiff University, UK, 2004. 
E. C. Titchmarsh. Eigenfunction expansions associated with second-order differential equations. Part I. Second 

Edition. Clarendon Press, Oxford, 1962. 



